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We develop a practical method of computing the stationary drift velocity V and the diffusion 
coefficient D of a particle (or a few particles) in a periodic system with arbitrary transition rates. We 
D ■ solve this problem both in a physically relevant continuous-time approach as well as for models with 

discrete-time kinetics, which are often used in computer simulations. We show that both approaches 
yield the same value of the drift, but the difference between the diffusion coefficients obtained in 
each of them equals \V 2 ■ Generalization to spaces of arbitrary dimension and several applications 
of the method are also presented. 

o 



I. INTRODUCTION 



CO i Investigation of diffusive transport is of highest importance in many areas of physics and related sciences. The most 
fundamental characteristics of diffusion is provided by two quantities - the diffusion coefficient D and the drift velocity 
V. The knowledge of the latter is especially important in studies of non-equilibrium phenomena, where usually V =/= 0. 
This includes, among others, diffusion in disordered media subject to an external field or the so called molecular 
' "pumps" and "motors" f8H10||, responsible for transport of various chemicals in biological cells. However, no universal 
and practical theoretical method of determining V and D for arbitrary systems has been developed, and although 
Q ■ many solutions for some particular physical models have been proposed, in more complicated cases one often has to 
O ' resort to approximations or numerical methods. 

The first, to our knowledge, attempt towards determining V and D in an arbitrary system was made by Derrida 
| jnj, who considered a one dimensional, periodic lattice of a period L and arbitrary hopping rates between nearest- 
. neighbor sites. He managed to give exact expressions for the velocity and the diffusion constant as functions of the 
hopping rates. He then extended his method |T^ ] to a periodic <i-dimensional system, but his solution turns out to 
be rather complicated. An alternative method of calculating V was also developed by Kehr et al. Q. Of the two 
quantities, V and D, the latter is of course much harder to find. A surprisingly simple formula for one-dimensional 
systems with transition rates satisfying the detailed balance condition was derived by Lyo and Richards fujfi , and 
recently independently reinvented by Kutner p4| and Wichmann p5[ , 




d=\t*Y,i^\ . (i) 



where L is the period of the lattice, Pj q denote the equilibrium site occupancy probabilities, and T~* are transition 



rates from site j to j + 1. However, this equation is valid only when the random walker hops between nearest- neighbor 
sites. A method of determining D in arbitrary periodic systems satisfying the detailed balance condition was then 
proposed by Braun and Sholl [n6|. Description of several other techniques, developed for some special types of random 
environments, can be found in review articles @-|§[L7|- 

All the above-mentioned methods are based on the following scheme: 



1. Reduce the infinite system to a single elementary cell with periodic boundary conditions. 

2. Write down the master equation. 

3. Using it calculate the steady-state properties of this system; in particular, find the steady-state site occupation 
probabilities 0,pUIi-|il| or propagators El . 
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4. Given these quantities, find V and D. 



The third step is here critical and can be carried out explicitly only for relatively small systems or for models possessing 
some special properties, e.g, one-dimensional lattices with jumps restricted to the nearest-neighbor sites. The main 
goal of this paper is to develop a simple and, at the same time, general method of calculating V and D explicitly 
without performing detailed analysis of the steady state. In our approach the third step reads 

3'. Find the matrix A(fc) representing the Fourier transform of the master equation. 

which is a rather straightforward operation. The technique we propose is valid for both equilibrium and non- 
equilibrium systems. It can also be quite easily implemented in computer-algebra or numerical programming lan- 
guages. 

Following the first step, our considerations will be restricted to periodic systems, such as crystals or molecular 
motors. A general problem of diffusion in an arbitrary disordered system remains an open question. One might be 
tempted to treated it by taking the limit L — > oo. This method can sometimes provide a hint about anomalous type of 
diffusion in a given aperiodic system (Tt]]. However, except for some simple models HJlfJ, it has not been established 
whether the limits t — > oo and L — > oo commute and thus, whether this approach is generally correct. 

We will pay a special attention to a one-dimensional system of finite period L with arbitrary jump rates between any 
of its sites whose distance is less than L. This is perhaps the simplest model of a periodic system where all elements 
of the transition rate matrix can take on nonzero values. Its investigation can be therefore carried out with relatively 
simple mathematical formalism. A solution to a general, multi-particle and multi-dimensional problem requires the 
same mathematical methods, but the notation would have to be complex. What is even more important, any periodic 
system with finite number of sites in its elementary cell can be mapped on such a one-dimensional system; geometry 
of the original problem is relevant mainly in constructing the matrix A(fc) and can be taken into account quite easily. 
Once the explicit form of A(fc) has been found, one can compute V and D using the methods we derive here for this 
simple one-dimensional system. 

In our calculations we will use both continuous-time and discrete-time formalisms. One could expect that since 
we calculate V and D in a stationary state, where all quantities are independent of time, both approaches should 
yield the same result. However, as was found by Derrida |]TT|| , the form of the diffusion coefficient as a function of 
the transition rates can depend on whether time flows continuously, or not. This is an important issue, since using 
discrete-time models belongs to favorite techniques of computer physics (in particular the cellular-automata p0[ , 
and the exact-enumeration || methods), and one needs to know whether results obtained in this way are correct. 
Derrida Jll| showed that in his (grossly simplified) model the relation between the diffusion constants obtained in 
the continuous-time (D) and discrete-time (-D D ) formalisms reads D — D B = \V 2 . We will show that this relation is 
general and will prove it for any periodic system with arbitrary transition rates. 

The paper is organized as follows. Using the Fourier transform of the master equation, see Refs. (]H p^7|j2l| l , in Section 
[n| we analyze a one-dimensional system with arbitrary period L and transition rates Tji . We find a simple technique 
of calculating V and D in terms of the 3 lowest coefficients of the characteristic polynomial of the Fourier-transformed 
transfer-matrix, A(fc). The most important properties of its spectrum at k = has been alre ady desc ribed in Refs. 



l|,[22f , but only for systems obeying the detailed balance condition. Therefore, in Sections II A and [I B we analyze in 



detail the spectr um of A(k) in a general case, assuming that the time is a continuous or a discrete quantity, respectively. 



Then, in Section II C we compare the two approaches and explain the different forms of diffusion coefficients derived 



in each of them. In Section |IH| we generalize our approach to systems in arbitrary space dimension, and in Section IV 
we present a particularly simple method of calculating V and D for one-dimensional systems with transitions between 
the nearest-neighbor sites. In Section [v] we derive explicit forms of V and D in some commonly used models. In 
particular, one of the examples studied there explains how to apply our technique to many-body problems. Finally, 
section [Vj is devoted to conclusions. 

II. GENERAL CASE OF A PERIODIC ID SYSTEM 

A. Continuous-time formalism 

Consider a one-dimensional lattice with its sites located at i„, n £ J\f. At time t = we put a particle at site 
xq = 0. The particle can then jump between the lattice sites. Transitions are assumed to represent a continuous (e.g., 
a Poisson) Markov process in time. The (constant in time) transition rate from a site x n to a site x m will be denoted 
by r mn . We assume that the system is periodic in space and denote its period by L > 1, 

^n,m T mn ^ m+L ,n-\-L j — -^m+L *^n+L- (2) 
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Of course transition rates between two different sites cannot be negative, T mn > if m ^ n. For simplicity we also 
restrict our considerations to the case where direct transitions between sites m and n are allowed only if \m — n\ < L. 

Let P(n, t) denote the probability density of finding the particle at site x n at time t. The evolution of this quantity 
is governed by the master equation 

dP{ ^ l) = FnrnPim, t) - T mn P{n, t)] , (3) 

m 

and the initial condition reads 

P(n,0) =8 n , . (4) 

We can now treat the whole lattice as if it consisted of L sub-lattices and define the probability densities Pi (j, t) of 
finding a particle at a given sublattice I at site n at time t as 

Pi(n,t) = P{n,t)5^, (5) 

where , is a generalized Kronecker's delta, 8^ tl = 1 if n = I (modL) and 5^ = otherwise. We can now rewrite 
(|J) as a system of L linear differential equations for Pi(n, t), I = 0, . . . , L — 1, 

J2 [r,, I+i n fJ (n + i,t)-r l+iiI fl(n,t)]. (6) 

j = l-L 

The major advantage of (|^) as compared with (Q) is that its coefficients do not depend on n. Therefore we can 
calculate its Fourier transform, 

E [r M+j e^-^p i+i (fc,t)-r i+jV PKfc,t)], (7) 

' ' j=l-L 

where Pi(k,t) = J2 n exp(-te n )P;(n, t). 

The system of equations (Q) can be written in a compact form using an L x L matrix Aij(k), 

dPl(k,t) -j^AimP^t), (8) 



dt 

3=0 



where 



Y lje iH*j-*i) + Tt J+L c ik( - x ^+ L - x '\ j < I 
Aij(k) = { Tije 1 ^"^ +r lj - L e lk ( x J- L - x '\ j > I (9) 

We will impose only one restriction on the form of transition rates Ty. We will demand that A(0) be irreducible (by 
a permutation of indices) 24 1, i.e., in the stationary state the particle can be found at any of the sublattices defined 
in (^) with a probability > 0. Physically this condition means that in the long-time limit the system does not split up 
into several non-interacting subsystems. Mathematically irreducibility means that A(0) has exactly one eigenvector 
whose all components are strictly positive. Our approach is thus quite general — we do not require that the transition 
rates should satisfy the detailed balance condition, the Fourier-transformed transition rate matrix, A(/c), need not be 
symmetric or even diagonalizable at k = 0, and some of its eigenvalues can be complex. 

A general solution to (|J) reads p5| 

L-l 

Pi(k,t) = J2T lj {k,t)exp[X j (k)t], 1 = 0,..., L-l, (10) 

3=0 

where the coefficients 2y(fe, t) are polynomials in t and can be determined using the initial condition. The degree of 
Tij(k,t) is smaller than the multiplicity of Xj(k). We assume that the eigenvalues Xj(k) are ordered in accordance 
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with the descending magnitude of their real parts at k = 0, i.e., j < I 
Yli=o t)> there is P(k, t) = J2i^=a Pi(k, t), and so ( |iT)| ) yields 



5R(Aj(0)) > ft(Aj(0)). Since P(n,t) 



L-l 



P(k,t) = J2hi(k,t)exp[X l (k)i\, 



(11) 



1=0 



where hi(k, t) = J^Lq 1 Tji(k, t) are polynomials in t (actually hi(k, t) can depend on time only if Xi(k) is degenerated). 
In a one-dimensional system the stationary drift velocity V and the diffusion coefficient D are given by 



V 



hmM 

t^oc t 



— lim i 

t — >oo 



dP(kA) 




dk 


k=0 


t 



(12) 



D = lim 

t—>oo 



(* 2 ) - {x? 
2t 



= lim 

t — >oc 



d 2 P(k.t) 




( dP(k,t) 


y 


dk 2 


+ 

k=0 


\ dk 


fc=0/ 



2t 



(13) 



where (f{x)) = J f(x)P(x,t) dx. Thus, if we could calculate the eigenvalues of A(fc), we would be able, in principle, 
to calculate V and D. The diagonalization of A(fc), however, is a formidable task feasible only in some special cases 
(e.g. for small L or when special relations have been imposed on its elements). Fortunately, as we will see below, 
if our aim is restricted only to determining V and D as functions of the transition rates Tij , we need not calculate 
explicitly even a single eigenvalue of the transition rate matrix! 

Let n be a constant such that Vj \i > \Ajj(0)\ and let Q denote an auxiliary matrix, Qjj = Aj;(0) + /!<5j-;. Since 
A(0) is irreducible, so is Q. Moreover, because all elements of Q are nonnegative, we can apply to it the Frobenius 
theorem |2^,|24j and conclude that Q has a positive eigenvalue q, which is a simple root of the characteristic equation, 
and the moduli of all other eigenvalues are at most q. Because the dominant eigenvalue of an irreducible, nonnegative 
matrix lies between the largest and smallest column sums ]23|p4| ], and in our case these sums are all equal to /i, we 
find that the dominant eigenvalue q = fi. Since the spectrum of Q is shifted, with respect to the spectrum of A(0), 
by (i, we conclude that for k = the matrix A(fc) has exactly one dominating eigenvalue Ao(0) = and the real parts 
of all other eigenvalues are negative, 



A o (0) = > R(Ai(Q)) >,...,> 3t(Ai_i(0)). 

Thus, in the limit t — > oo this single eigenvalue dominates the sum in r.h.s of ([i 
of V and D we can therefore employ a simple approximation 

P(k,t) = h (k,t)exp[X (k)t]. 

Moreover, since J2 n P{ n ^) — 1> there is 

V t > P(0,t) = 1, 

which, owing to (|l4|), implies 

h o (0,t) - 1. 

Upon inserting ( |l5| ) into (111), ( |l3| ) and using ([14]), dljj ) we conclude that 



V = i 



D = - 



k=Q 



dk 
ld 2 X 



2 dk 2 



(14) 

at k « 0. In determining the forms 

(15) 

(16) 
(17) 

(18) 
(19) 



k=0 



Let W(x) denote the characteristic polynomial of the matrix A(fc). Let Cj(fe) denote its coefficients at x- 7 , j 
0, . . . , L. We thus have 



4 



V fc VF(A o (fc))=^ Cj (fc)[A o (fc)P=0. 

3=0 

On differentiating it with respect to k and using the first part of we arrive at 



(20) 



dk 



k=0 



(21) 



where we used a short-hand notation Cj = Cj(0) and c'j = dcj /dk\k=o- Note that c\ ^ 0, which follows from ([j^ 
Differentiating (pQ) twice yields 



3 2 A 



dk 2 



4 + 2c 2 (X' ) 2 + 2c[X' 



k=o c i 

where Ag = dXo(k) /dk\k =0 and c ' = 9 2 c /3fc 2 |fc =0 . We thus finally arrive at our major result 



v 



D 



H) 
'■ — ; 
Cl 



c ' - 2c 2 V 2 - 2ic[V 
2^ ' 



(22) 

(23) 
(24) 



Functions Cj{k) can depend on k and the transition rates in a very complicated way. We were able to find only two 
general properties, both following immediately from jl^): cn = and c\ ^ 0. Explicit forms of c , c ', ci, c[ and C2 
for some particular models will be given below, in sections IV and M. 



B. Discrete-time formalism 



The discrete-time formulation of the problem is basically similar to the continuous one, but there are some major 
differences, too. For simplicity we will use the same notation as in the previous section, but one should remember 
that almost all functions employed in the discrete-time formulation of the problem differ from those we dealt with 
in Section II A. Where it will be necessary to compare quantities computed within each approach, we will attach a 
superscript "D" to the quantity derived in the discrete-time formalism. 

In the discrete-time version of the problem the master equation (expressed in terms of the probabilities Pi(n,t) of 
finding a particle at time t at a site x n belonging to a sublattice /) reads 



L-1 



Pl(n,t+l)=P l (n,t)+ [ T i,i+i p i+j( n +.h t )- T i+j.i p i( n ^)] 

■j = l-L 



(25) 



where I = 0, . . . , L — 1, and are dimensionless probabilities satisfying the usual condition Vy < Tji < 1. Note that 
in the continuous formulation of the problem Tji were unbounded from above, dimensional quantities (of dimension 
[T -1 ]) and we called them "transition rates" . Last, but not least, in the present approach time t assumes only integer 
values. 



Upon taking the Fourier transform of (25) we arrive at 



L-1 

P l (k,t + i) = P l (k,t)+ J2 \?i,i+^ lk(xi+3 ~ Xl) Pi+ 3 (M) - ?i+i,iPi % *) 

j=l-L 



(26) 



which can be rewritten using a stochastic matrix AR(fc) 



L-1 



3=0 



(27) 



where the only difference between A D and its continuous-time counterpart A lies in their diagonal elements 

Af j (k) = A lj (k) + 5 lj . (28) 
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Therefore the eigenvalues Xf of A D are related to the eigenvalues A; of A through a simple formula 



Af(fc) = Aj(fc) +1, 1 = 0, 



,L-1. 



Applying the theorem of Frobenius |23| , P4|] to A D (0), which is a real, nonnegative matrix, we conclude that 

A? (0) = 1 



(29) 



(30) 



and the moduli of all other eigenvalues do not exceed 1. However, in contrast to the continuous-time formalism, now 
there can be s > 1 eigenvalues, say Xf m (k), m = 0, . . . , s — 1, such that \Xf m (0)| = 1; for convenience we assume that 
Aj^ = Xq . Although this can happen only if all diagonal elements of A D (0) vanish we do not exclude this 

exceptional case from our considerations (see Ref. pl|). The Frobenius theorem ensures us that all these dominating 
eigenvalues are distinct, and so in the normal Jordan representation of A D (fc) the size of the corresponding Jordan 
blocks is 1. This suffices to assert that for k « and t — > oo 



P(k,t) 



s-l 
m=Q 



(fc)[A°(*0f 



(31) 



where h m are some functions of k. 

Because equation ( [l6|) is valid both in continuous-time and discrete-time formalisms, upon comparing it with (|3l| ) 
we conclude that except for ho(k) all other coefficients h m (k) in (|3l]) must vanish at k = 0, 



M0) = 1, Mo)=o, m = l, 



1. 



(32) 



Actually, since t is an integer, this is not a trivial state ment; a proof is based on the fact that all dominant eigenvalues 
AS are distinct roots of the equation x s = 1, see Refs. 
(pp|) and (|3^) we arrive at 



jjjpil- Upon inserting (^) into (|l^) and (|l^) and then using 



V L 



D 



dX»(k) 



Ok 



k=0 

dk 2 



( d\%(k) 
\ dk 



*(0,i), 



fc=0 



where 



*(*:,*) 



dh m ( dX 



m—l 



dk \ dk 



9>l^ 
dk 



A m (A m ) 



D \t-l 



(33) 
(34) 

(35) 



Note that, by definition, ^(k,t) = if s = 1, i.e., if Ap (0) is the only dominating eigenvalue of the transition matrix 
A(0). Moreover, because the existence of the diffusion constant is guaranteed by the central limit theorem, \l/(0,i) 
must be independent of time, at least in the limit t — » oo, for any s. This, in turn, requires that even if s > 1 



f(0,t) = 



(36) 



(see the comment under eq. (|32|} ) . Therefore, to determine V and D one needs only to investigate the properties of a 
single eigenvalue A^(0) at k w 0, i.e., eq. (J3l|) can be replaced with 



P(k,t) = h (k)[X»(k)] t 



(37) 



Finally, using ([lj), (|l9|), (p9|), and ( |33"| ) - (^) we conclude that the drift velocities and diffusion coefficients obtained 
in continuous and discrete formulations of the problem are related to each other by simple, general formulae 



D D 



V, 
D 



±V 2 

2 ' ■ 



(38) 
(39) 



Equation (|39| ) holds, however, only if the interval r between successive jumps in the discrete formulation of the 
problem equals 1, which has been assumed in our calculations for simplicity. For a general value of r the relation 
between D and D u reads 



D u = D- ±tV 



(40) 



which has a correct dimensional form (all its terms have a dimension [L 2 T 



G 



C. Comparison of the two formalisms 



Comparing equations (|15|) and ( |37| ) we find that the main mathematical differences between the continuous- and 
discrete-time formalisms are related to different functional forms of P(k,t) at k s=s 0. While in the continuous-time 
approach this function is given by an exponential, in the discrete-time model we deal with a power function. In 
particular, notice that d 2 X t (k)/dk 2 yields a term proportional to t(t — 1) = t 2 — t. This gives rise to an additional 
term linear in t, which is responsible for the difference between ( |34| ) and (|l9|), and hence for the term \tV 2 in (|40|). 

The presence of the time unit r in ([lO]) suggests that the key to understanding the physical reasons for the difference 
between D and D° lies in the dimensional analysis. In the continuous-time approach the transition rates Tji are 
dimensional quantities that scale with time as [T -1 ]. Therefore, multiplying them all by some positive constant a/1 
corresponds to changing the time unit, and so the resulting diffusion coefficient D ncw will be equal to the product 
of a and the original value of the diffusion coefficient, Z? id- Similarly, V ncw = ctVoid- This explains why in the 
continuous-time approach we could safely assume r — 1. However, in the discrete formalism there is no such a simple 
relation between the time interval r and the jump probabilities Tji, which are dimensionless. Suppose, for example, 
that at times t = 0,r, 2r, ... we toss a coin. Changing the frequency of tossing, or 1/r, will have no impact on the 
probability of the coin falling heads up. The proper way of taking the continuous-time limit in the discrete-time 
model is to assume that the jump probabilities T^j are related to the continuous-time transition rates Tji by a simple 
formula 1^ = rTji and taking the limit r — » 0. In other words, to get the continuous-time limit, we need to apply the 
"alpha-transformation" to the transition probabilities with an infinitesimally small value of a. Equation (^) shows 
that D D is actually a sum of two terms — one that scales linearly under the "a-transformation" (D), and the other 
one which scales quadratically (—^tV 2 ). In the limit r — > the continuous-time diffusion coefficient D is thus of 
order 0(t), while —\tV 2 is of order 0(t 2 ) and can be neglected. Consequently, as could be expected, the relative 
difference between diffusion coefficients calculated within each approach vanishes, (D — D°)/D — > as r — > 0. 

A decrease of the diffusion coefficient in the discrete-time formalism can be also interpreted as a consequence of 
the fact that a discrete process tends to be "less random" that a continuous one. This is clearly seen in a limiting 
case of a one-dimensional lattice (with lattice constant a and time unit r) where all probabilities of jumping to the 
left vanish (r*~ = 0) and all probabilities of jumping to the right are equal 1 (rj* = 1). A discrete process with this 
choice of probabilities is completely deterministic. At each time interval the diffusing particle with probability 1 hops 
to the right, hence the diffusion constant vanishes. On the other hand, if we consider a continuous-time process with 
= and = 1, wc can see that the motion of the particle is now by no means deterministic. Although we know 
that on average there will be one jump to the right per unit time, we do not know when it will actually occur. This 
uncertainty introduces randomness to the process, and the corresponding diffusion coefficient equals a 2 /2r. 

Note finally that because the left-hand side of (|4^) represents a diffusion coefficient, it cannot be negative, and so 

D > \tV 2 . (41) 

This has several interesting consequences. First, if in a periodic system there is a stationary drift (V / 0), then the 
continuous-time diffusion coefficient D is bounded from below. That such a bound exists can indirectly imply bounds 
on other physical quantities, e.g., the maximal force exerted by a molecular motor p6| . Second, suppose that (|4C|), 
and hence (f4l|), is also valid in infinite (aperiodic) random systems. Whenever the diffusion is sublinear, i.e, whenever 



lim t _ >00 ((x ) — (x) )/t = 0, there is D = 0. Equation (40) would then imply that V = 0, or lim^ 00 (x)/i = 0. This 



would mean that whenever diffusion is sublinear ("subdiffusion'), the drift is also sublinear ("subdrift") or vanishes 
altogether. The asymmetric hopping model with bond disorder can serve as an example of an infinite system where 
such relation is actually observed M. 



III. DIFFUSION IN ARBITRARY SPACE DIMENSION d 

Suppose we want to calculate the drift velocity V and the diffusion tensor D in a system with d Euclidean co- 
ordinates. If the system can be divided into a finite number of subsystems with constant transition rates between 
each two of them, our major results (14) - (0) and Q37]) remain valid irrespective of the geometry of the system. 



Such division is always possible for periodic systems with finite number of sites in the elementary cell. The matrix 
A(fc) is the most "sensitive" to the geometry of the system under consideration. In calculating its explicit form one 
can use an equation similar to ([)]), remembering, however, that k and x n are now vectors in a d-dimensional space. 
One should also take into account all possible transitions between sublattices, and this can be done by considering 
all possible jumps starting at any site belonging to some elementary cell and ending in the same cell or in one of its 
nearest-neighbor cells. 
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The components of the velocity vector V and the diffusion tensor D are given by 



V„ 



lim ^ 

t^oo t 



.1 dP{k,t) 

urn i — 



t^oo t dk u 



D 



lim 

t — >oo 



2t 



fc=0 

lim — 

t— >oo 2t 



d 2 P(k,t) dP(k,t) dP(k,t) 



dk^dk a 



dk,. 



dk n 



fc=0 



Using (H) - @ and @ we conclude that 



V u = VP = 



dk u 



fe=0 



D, 



D 



1 d 2 X 



2 dk^dku 



k=0 



dk^dkv 



dkfj, dk a 



(42) 
(43) 

(44) 
(45) 
(46) 



fc=0 



where [j,,cr = 1, . . . , d. Actually, since log P(k, t) rj \o(k)t is a generating function for cumulants of P(x, t) (see p2f), 
the first two of the above formulae are a natural consequence of (15) and (|l7|). 

Just as in the one-dimensional case any component of V or D can be expressed in terms of derivatives of the three 
lowest coefficients of the characteristic polynomial of the matrix A(/c) at k = 0. In particular, 



D 



2c 2 V^ CT - * (c^Va + 4' 



(J,<7 



2 Cl 



fc=0 



D 



D, 



1 



2 

where we used a short-hand notation = df /dk^ 



(47) 

(48) 
(49) 



IV. THE NEAREST-NEIGHBOR TRANSITION RATES 



If in a one-dimensional system only transitions to the nearest-neighbors are allowed, and if the distance between 
nearest- neighbor sites is equal 1, our formulae (E3h and (|24|) can be simplified by noting that in this particular case 



(50) 



c' Q = iL\ 


(L-1 


L-1 

n 




\3=0 


4=0 






L-1 


c» = L 2 




4=0 



(51) 



ci = 0, (52) 

where Tj* and are the transition rates to the right and left from site j, respectively; r~* = ^j+i.j and = I^-x^-. 
Thus, to determine V and D we only need to compute the two terms of the characteristic polynomial of A(0): c\ and 

C2- 

Although ci and ci could be, at least in principle, found by calculating det(a;/ — A(0)), / being the identity matrix, 
this could lead to serious computational overhead even for L of order 10. A more efficient method exploits the fact 
that ci are polynomials in and Tj* , j = 0, . . . , L — 1, and can be expressed as 
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L-1 

C l= E U ^^^n^Mll^ASj}) (53) 

where Z 6 {1, 2}, j m £ {0, 1}, S n € {0, 1}, and ^>;({7j}, {<5j}) = if at least one of the following conditions is satisfied: 

E m =o(7m +6 m )^L-l, 

3 m 'Jrn — : — 1? 

3 m 7m = <Wi = 1; (54) 
otherwise "0({7i}j = 1- For example, for L = 3 there is 

c i = r r r r + r r r tT + rrrr 

+1T17 + iriy + iYiy 

+r^rr + r-rr + r-rr, 
C2 = r^ + rr + rr + r^ + rr + r^. (55) 

Although formally the sum in ( p3| ) consists of 2 2L terms, in practice only for a few of them ^({7.?}) 0- in 

particular, for a given value of L the sum in ( |53] ) consist of L 2 non- vanishing terms for c\ and (L 4 — L 2 )/12 terms for 
C2- Thus it should be possible to calculate D algebraically for L of order 20, and numerically for L of order at least 
100. It is worth noting that since cj can be computed as sums of positive values, by using ( p3| ) rather than calculating 
the determinant of xl — A(0) we can avoid large numerical errors which often appear when computing determinants 
of large matrices. 

It is not difficult to show that our general formulae (|2^) and (^J) are consistent with eq. ([!]) derived for a general case 
of ID systems at equilibrium with nearest-neighbor transitions. To this end it suffices to notice that the equilibrium 
site occupation probabilities P° q can be expressed in terms of the jump probabilities as 



' • 1 .- 1 ..• 1 ,< ' j i 1 .- 1 / I--- 1 .,. /. • : 



(56) 



On inserting it into (|l|) and then taking into account (53) and the trivial condition V = one arrives at (p4|). 

Another interesting consequence of ( p0|) - (|5^) is that diffusion in systems with transitions to nearest-neighbours is 
always normal. To see it note that since c\ ^ 0, D is always finite and thus no superdiffusion is possible. The other 



type of anomaly, subdiffusion, would require that D = 0. Owing to (41) this would imply V — which, upon taking 
into account (pi|), would require Cq = 0. However, the explicit form of Cq, see (|5l|), guarantees that Cq > except for 
a degenerated case where the diffusing particle is confined to a finite region of the lattice (TJ" = Tp = for some j 
and I). 

One final remark. Some authors p[ JlO|Jl4] |, when considering a one-dimensional system with nearest- neighbor tran- 
sitions, prefer to reduce the problem to diffusion on a finite ring and investigate the probability current J rather than 
the drift velocity V. The former is defined as J = PjTj* — Pj+iT^_ 1 and in the steady state does not depend on 
the site j it is measured at. For single-particle systems J and V are related to each other through a simple formula 
J = V/L [jflj. In this case our approach can thus be applied for calculating both V and J. 

V. APPLICATIONS 

A. The case L = 1 and L = 2 

Applying our approach to a particle diffusing in a one-dimensional system with a lattice constant a = 1, time unit 
t = 1 , and the period L = 1 we immediately arrive at a well known result 

V = T^-T^, (57) 

# = (r^ + rn/2, (58) 
d d = [r^ + - (r^ - r^) 2 ]/2. (59) 
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For L = 2 we find 

V = 2 rrry ~ rrIT , (60) 

D = 2 Errr+rrrr_^ (61) 

pn _ 8 my+rrrr .!+g^ i (62) 

where S = T^* + + + Ttf. Note also that the solution for L = 3 can be easily derived using (|55|). 



B. The sawtooth potential of arbitrary period in an external field 

Consider a one-dimensional system of period L with site energies E n L + j = je, where n is an integer, j = 0, . . . , L— 1, 
and e > is a constant. Such a pattern is know as a discrete sawtooth potential ||[|- We assume that the transition 
rate from a site j to j + 1 is given by 

iy = 6exp(-/fe/2), (63) 

and the rate of jumping from j to j — 1 reads 

/ 6- 1 cxp(/3e/2), j ^ (mod L) , . 

3 \b- 1 exp(/3e{l/2- L)), j = (mod i) l ° J 

where /3 = l/k^T is the Boltzmann factor and b = exp(/3F/2) represents a bias due to an external force F. 
Note that with this choice of the transition rates, for b = 1 the system satisfies the detailed balance condition 
T^ +1 exp(-l3E J+1 ) = iy exp(-/3^). 

Using ( p3[ ) and some combinatorics one can prove that for L > 2 

c' Q =iLR L (b L -b~ L ) (65a) 

4 = L 2 R L (b L +b- L ) (65b) 
ci = G- L - 1 5 L (G) + fl 2L G L_1 Si,-i(G _1 ) (65c) 
c; = (65d) 



c 2 = G^^Z L _ 2 (G) + R ZL G^Z L ^{G-') (65e) 

where i? = exp(— j3e/2) controls the anisotropy of the potential, G = S m (x) = J 3 ' 2 ' 7 = [x 2m+2 (1+m — mx 2 ) — 

x 2 ]/(x 2 — l) 2 , and Z m (x) = | X)Jlo( m + 1 ~ j)0 + 1) (j + 2)x 2: '. Actually we were able to prove only (65a)- 
and (p5a) is a conjecture based on the form of c 2 derived for L = 2, . . . , 20 



The stationary drift velocity V and the diffusion constants D and D D can be now calculated using (J2^), (^J), and 
( pp| ) . The resulting velocity has already been studied in Ref. [0 . The diffusion coefficient D calculated for various 
values of R and L as a function of the bias b is depicted in Fig. [1. For R = 1 there is e = 0, and so for any j we have 
Tj^ = b and iy = b~ l . Consequently, the effective period equals 1. Using fl5S| ) we conclude that D = ^(fr + fr -1 ). For 
R < 1 the behaviour of D becomes more complicated. For a large bias to the right, b 3> 1, the jumps to the left are 
so rare that practically they become irrelevant. We thus have = bR, w 0, and so D w 6i?/2 irrespective of L. 
For a strong bias to the left, b <C 1, the jumps to the right can be neglected, but the particle has to jump over a large 
potential barrier located at sites j = . . . , —L, 0, L, . . . whose height is proportional to L — |, see (f34|). Therefore, in 
this regime D is a quickly decreasing function of L. These two limiting solutions match in the intermediate regime, 
which can be roughly described as 1 < b < 1 / R. 



C. Two particles in a sawtooth potential on a ring of length L = 4 



Consider two particles diffusing in a sawtooth potential on a ring of length L — 4. We assume that any site j can 
be occupied by only one of them (the hard-core interaction) . For simplicity we also assume that the distance between 
the particles cannot exceed L — 1, so that we can reduce the system to a ring consisting of L sites. 

Each state of the system can be described as a pair of integers, (n, m), where n, m = 0, . . . , L— 1 denote the currently 
occupied sites. As n ^ m, there are L(L — l)/2 = 6 different states. Our two-particle system is thus equivalent to a 
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6-state system with one "virtual" random walker. These states are, in order, (0, 1), (0,2), (0,3), (1,2), (1,3), (2,3). 
The distance between a state (ii,ji) and («2,j2) equals ii — i\ + ji — j\. For example, the distance from (0,1) to 
(0,2) is 1. The transition rate between (ii,ji) and (12,32) vanishes unless \%2 — ii\ + |j2 — Ji = 1, and in this case 
a transition from («i, ji) to (12,^2) corresponds to a single jump of one of two diffusing particles. For example, the 
transition rate from (0,1) to (0,2) is T^*, and from (0,1) to (0,3) is 0. Consequently, if the lattice constant a = 1, 
the matrix A(k) takes on the following form, 



A(fc) 





ITe ifc 








rr e - Ife 





rre~ sfc -rr - 


- rr — r — r 


rre 8fe 







rre- ife 





iye- ifc 


-rr - rr 





rre 8fe 








IYe- ifc 





-rr - rr 


rre 8fe 










-ik 

1 e 


rre-* fe 


-rr - rr - rr - rr 


rre Ife 





ITe ifc 








rre" 8fe 


-rr - rr 



(66) 



The drift velocity V and the diffusion coefficient D can be found for arbitrary transition rates, but the results are 
very lengthy. However, in a particular case of the sawtooth potential the transition rates Tj* and T*~ are given by 
relatively simple formulae ( |63| ) and (|64|), respectively. Then 



c o 
4 

Cl 



C2 



- 6 (G 8 



8iG 
32G 
G~ 5 (G 2 + 
+R s (9G e - 
+R 16 (G 4 - 
8iG- 5 (G 8 
G~ 4 [37G 8 
+i? 8 (23G 6 
+i? 16 (3G 4 



2G 2 - 
-2G 2 



-i? 8 )(^ 
+ R 8 )(R 
1)[12G 8 + 15G 6 
I- 13G 4 + 17G 2 4 
- 2G 2 + 5)] 
-i? 8 )(i? 8 + 4G 2 
+ 69G 6 + 58G 4 - 
+ 49G 4 + 53G 2 
+ 8G 2 + 9)], 



1)(G 2 ^ 
-1)(G 2 
f 10G 4 + 
5) 

+ 3) 

- 31G 2 + 
f 19) 



-1) 
M) 
5G 2 



(67) 



where G and R were defined below eq. (65). 

Just like in the previous example, V, D, and D° can be now determined using (|23|), (p4]), and (^). The properties 
of the velocity V will be studied in detail elsewhere (see Ref. pTjl ). Here, in Fig. g we present the ratio of the diffusion 
coefficient calculated for a 2-particle system (D 2 ) to that calculated for a single-particle system (Di), for L = 4, as a 
function of the bias b. This ratio measures the change of diffusivity due to interactions with the second particle. As 
could be expected, D2 < D±, i.e. the presence of the second particle decreases the diffusivity of the first one. This 
effect is not strong, however, for we find that D2 > \D\- As in the previous example, for large bias b, the ratio D2/D1 
becomes independent of the potential anisotropy R. 



VI. CONCLUSIONS 

We have developed an effective technique of calculating the drift velocity V and the diffusion tensor D in arbitrary 
periodic systems. A novel feature of our approach, as compared to previously proposed methods p|-|l6||, consists in 
the fact that we need not examine in detail the steady-state properties of the system. In particular, we do not need 
to solve linear equations of relatively high order to find the stationary probability distribution of the random walker 
over different sublattices. Instead we calculate V and D directly. The only quantity we really need in order to carry 
out our calculations is the matrix A(fc). Its explicit form, however, can be found almost trivially; all one needs to 
know are transition rates and distances from any site of an elementary cell to any other site in the same or adjacent 
cell. Once the form of A(fc) has been determined, one should calculate the three lowest terms of its characteristic 
polynomial. Then V and D can be determined using our formulae ( |47| ) and (Q), respectively. 

Our approach is very general and can be applied practically to any periodic system in an arbitrary space dimension. 
It is particularly well suited for calculations employing computer-algebra systems, e.g., Maple or Reduce, or for 
numerical analysis. In this context one should also note that the coefficients q of the characteristic polynomial of 
A(0) can be always expressed as polynomials in Aj;(0), j / /, with positive coefficients, and thus can be calculated 
numerically with extremely high accuracy. While this is not necessarily true in the case of their derivatives, we believe 
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that our approach enables one to calculate V and D even for systems where the corresponding steady-state problem 
is numerically ill-conditioned. 

We were able to prove a general relation between the diffusion constants calculated using continuous-time and 
discrete-time approaches. Our analysis shows that in any periodic system (including, for example, lattice gases 
studied in Refs. they differ by a term \V 2 , where V denotes the steady-state drift velocity along a given 

direction. Since the limits t — > oo and L — > oo were shown to commute, at least in some simple models |2|Jl9(|, we 
expect that this relation holds also for infinite random systems. It would be interesting to know whether such a simple 
and apparently universal relation can be applied in contexts other than those studied here. 

We also showed how our technique can be applied to a simple many-body problem. While, owing to mathemat- 
ical complexity, we do not expect that in this way one will be able to calculate explicitly transport coefficients in 
non-equilibrium systems containing more than a few particles, our conclusion about the universal difference between 
diffusion coefficients calculated in continuous- and discrete-time models should still hold. This conjecture is based on 
the fact that diffusion of several particles can be interpreted as diffusion of a single particle in the corresponding mul- 
tidimensional phase space, and for the latter problem the relation between D and D° has been established rigorously 



in Section HI. Thus, for the future work, two problems are of primary importance: the role of periodic boundary 
conditions in the limit of asymptotically infinite period, L — > oo, and detailed analysis of possible applications of our 
method to many-body systems. 
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FIG. 1. Diffusion coefficient D of a single particle in a sawtooth potential (|63|) and ( p4| ) for various potential periods L and 
anisotropy parameters 7? as a function of the bias b. Arbitrary units. 
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FIG. 2. The ratio of the diffusion coefficients D2 and Di calculated for a 4-site ring containing 2 and I particles, respectively, 
as a function of the bias b, for various values of the anisotropy parameter R. Arbitrary units. 
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